#include "Splines.h"

# define PI acos(-1)

double distance(double t1,double t2)
{
    double x1,x2,y1,y2;
    x1=sqrt(3)*cos(t1);
    x2=sqrt(3)*cos(t2);
    y1=2.0/3*(sqrt(3)*sin(t1)+sqrt(sqrt(3)*fabs(cos(t1))));
    y2=2.0/3*(sqrt(3)*sin(t2)+sqrt(sqrt(3)*fabs(cos(t2))));
    return sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
}

double curvelength(double t1,double t2)
{
    double tm=(t1+t2)/2;
    if(distance(t1,tm)+distance(tm,t2)<=1.000001*distance(t1,t2)) return distance(t1,tm)+distance(tm,t2);
    else return curvelength(t1,tm)+curvelength(tm,t2);
}

double Newton(double s,double l,double r)
{
    double m=(l+r)/2;
    if(curvelength(-PI/2,l)==s) return l;
    else if(curvelength(-PI/2,r)==s) return r;
    else if(r-l<0.0000001) return m; 
    else if((curvelength(-PI/2,l)-s)*(curvelength(-PI/2,m)-s)<=0) return Newton(s,l,m);
    else return Newton(s,m,r);
}

int main(int argc, char * argv[])
{
    double totallength=curvelength(-PI/2,3*PI/2);
    int N=atoi(argv[1]);
    vector<double> t,x,y;
    for(int i=0;i<=N;i++)
    {
        double tt=Newton(totallength*i/N,-PI/2,3*PI/2);
        double xx=sqrt(3)*cos(tt);
        double yy=2.0/3*(sqrt(3)*sin(tt)+sqrt(sqrt(3)*fabs(cos(tt))));
        t.push_back(tt);
        x.push_back(xx);
        y.push_back(yy);
    }
    Splines sx,sy;
    sx=interpolation(t,x,0,0,10,10,"ppForm","cubic","ncs");
    sy=interpolation(t,y,0,0,10,-10,"ppForm","cubic","ncs");
    fstream f;
    f.open("dataE.txt", ios::out);
    f << t.size() << endl;
    for (int i = 0; i < t.size(); i++)
        f << fixed << setprecision(6) <<t[i] << " ";
    f << endl;
    for (int i = 0; i < sx.coe.size(); i++)
    {
        for (int j = 0; j < sx.coe[i].size(); j++)
            f << fixed << setprecision(6) <<sx.coe[i][j] << " ";
        f << endl;
        for (int j = 0; j < sy.coe[i].size(); j++)
            f << fixed << setprecision(6) <<sy.coe[i][j] << " ";
        f << endl;
    }
    f.close();
}